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The dominant reaction pathway (DRP) is a rigorous framework to microscopically compute the most probable 
trajectories, in non-equilibrium transitions. In the low-temperature regime, such dominant pathways encode the 
information about the reaction mechanism and can be used to estimate non-equilibrium averages of arbitrary 
observables. On the other hand, at sufficiently high temperatures, the stochastic fluctuations around the dominant 
paths become important and have to be taken into account. In this work, we develop a technique to systematically 
include the effects of such stochastic fluctuations, to order ksT . This method is used to compute the probability 
for a transition to take place through a specific reaction channel and to evaluate the reaction rate. 



I. INTRODUCTION 



The theoretical investigation of the kinetics of rare conformational reactions of macromolecules represents a fundamental, yet 
very challenging task. In view of the large computational cost of performing Molecular Dynamics (MD) simulations of the long- 
time dynamics of molecular systems, alternative methods have been developed, which allow to sample directly the space of the 
reactive trajectories in large conformational spaces, without investing time in simulating thermal oscillations in the (meta)-stable 
statesCHl. 

In particular, the Dominant Reaction Pathways (DRP) approach |5 -8| yields by construction the set of most statistically 
significant reactive trajectories in the over-damped limit of Langevin dynamics. In such an approach, the reactive paths are not 
calculated by integrating the equation of motion of the system. Instead, they are obtained by minimizing a target functional, 
which is rigorously derived starting from the original Langevin equation. 

The main advantage of the DRP approach is that it allows to remove the time as the independent variable. Instead, the 
dominant path is calculated as a function of a curvilinear abscissa / which measures the distance covered in configuration space. 
This way, the problem of the decoupling of the time scales in the internal dynamics of molecular systems is rigorously bypassed 
and typically (9(100) path discretization steps are sufficient to characterize an entire conformational transition. In the same 
formalism, the time-dependent dominant pathway x{t) can be rigorously calculated a posteriori, from the trajectory x{l). 

A potential limitation of the DRP approach is that it emphasizes the role played the most probable reactive trajectories. These 
are smooth paths, defined as the maxima of the functional probability density ^[x{l)] for the system to make a transition from 
given reactant to given product configurations. Clearly, the real physical transitions never occur through such smooth dominant 
paths, but only through non-differentiable stochastic trajectories. However, it is possible to show that in the low-temperature 
limit the path probability density ^[x] is peaked in the regions of the functional path space surrounding the dominant reaction 
pathways. In such a temperature regime, each dominant reaction pathway can be considered as representative of a different 
reaction channel — see Fig. [T] — . Hence, if the reaction can occur through n different channels — i.e. there are n distinct 
dominant reaction pathways — then the probability of making a transition through the /— th channel can be estimated from the 
equation 

Prob. (i-i\\ reaction channel ) — ^ ^ , (1) 

where ^[xk{l)] is the probability density of the ^— th dominant reaction pathway, Xk{l). In this formula, the stochastic fluctuations 
around the dominant paths are completely neglected. 

The amplitude of the stochastic fluctuations grows with the temperature of the heat-bath — see FigJT] — . In general, we 
expect that in many practical applications it is necessary to go beyond the leading-order approximation ([TfTi.e. to account also 
for the effect of small stochastic fluctuations around the dominant paths through an expansion in the thermal energy ksT . As we 
shall see, to the next order in ksT the probability of a reaction channel can be cast in the form 
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FIG. 1 : Stochastic fluctuations around two different dominant reaction pathways in a two-dimensional transition. As the temperature raises, 
the stochastic fluctuations become large and eventually the two paths become statistically indistinguishable. 



where f[xk] is the o{kBT) correction to the probability of the reaction channel defined by the ^-th dominant path. 

The purpose of the present paper is to study the role of such stochastic fluctuations in the reaction kinetics. In the first part of 
the work, we provide a rigorous and practical method to compute the f[xk] coefficients, i.e. the contribution to the normalized 
path probability T[x] arising from the path integral over the small stochastic fluctuations around the most probable paths, to 
order ksT . In general, numerically performing such a calculation can be quite computationally demanding. However, we shall 
see that using a recently developed formulation of the Langevin dynamics at low-time resolution power 11411151 it is possible to 
drastically reduce its computational cost. 

In the second part of this work, we derive an expression which relates the reaction rate to the probability densities in config- 
uration space, evaluated in the DRP approach to order ksT . In order, to illustrate and test our method with a high numerical 
accuracy, in this work we choose to focus on thermally activated transitions in simple low-dimensional toy systems and leave 
the analysis of much more complicated molecular reactions to our future work. We shall show that, once thermal fluctuations 
around the dominant paths are correctly taken into account, it is possible to obtain accurate estimates for the reaction rates and 
equilibrium population ratios. 

The paper is organized as follows. In section |ll] and |lll| we review the DRP formalism. In the subsequent section |IV] we show 
how to compute the contribution to the reaction path probability arising form small stochastic fluctuations around the dominant 
paths. In sectionjVj we show that the efficiency of such a calculation can be greatly improved by adopting the low-time resolution 
formulation of the Langevin dynamics, which was recently developed in [14, 15 1. In section VI we test our calculation of the 



path probability in an analytically solvable model and we find that it gives very accurate results. In section VII we derive an 



expression for the reaction rates in terms of the path probability density calculated in the DRP approach, while in section VIII 
we test this formula by computing the reaction rate in some test systems. Our results and conclusions are summarized in section 

na 



II. PATH INTEGRAL REPRESENTATION OF THE OVER-DAMPED LANGEVIN DYNAMICS 

Let us consider a system defined on a generic J-dimensional configuration space. We shall assume that the dynamics is defined 
by the over-damped Langevin equation 

i=-pZ)Vf/(x)+r|(0, (3) 

where x is a point in the configuration space, P = l/^^r, Z) is the diffusion coefficient, U {x) is the potential energy function and 
\\{t) is delta-correlated Gaussian noise, satisfying the fluctuation-dissipation relationship: 

{r\\t')r\\t)) = 25'^Z)5(r -?') = \,d). (4) 

Note that in the original Langevin Eq. there is a mass term mx. However, for macro-molecular systems this term is damped at a 
time scale 10~^^ which much smaller than the time scale associated to local conformational changes. 
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The stochastic differential Eq. ^ generates a probabiHty distribution P{x^t) which obeys the well-known Smoluchowski Eq.: 

^P{x,t)=Dy[yp{x,t)^^yU{x)P{x,t)]. (5) 
ot 

Such an Eq. is often written in the form of a continuity equation, 

|p(x,0 = -V/(x,0, (6) 

where 

J{x,t) = -Z)[V + pVt/(x)] P{x,t) (7) 

is the so-called probability current. 

By performing the formal substitution P{x^t) = ^(x,^), the Smoluchowski Eq. ([s]) can be recast in the form of an 

imaginary time Schrodinger Eq.: 

-^^{x,t)=Heff'V{x,t), (8) 

where 

Heff = -DV^ + Veff{x), (9) 

is an effective "quantum" Hamiltonian operator and 

yeff{x) = ^ [iyU{x) f - ^ V2C/(X)) . (10) 

is called the effective potential. 

The conditional probability P{xf^t\xi) to find the system at the configuration Xf at time t, provided it was prepared in the 
configuration Xi at time ^ = is the Green's function of the Smoluchowski Eq., i.e. 

jP{xf,t\xu^)-Dy[yP{xf,t\xu^)^^yU{xf)^^^ (11) 

Formally, the conditional probability P{xf^t\xi) can be related to the imaginary time propagator of the effective "quantum" 
Hamiltonian ([9]): 

P{xf,t\xi) = ^-i(^(^/)-^(^'-))/r(x/,r|xO= e-^2iu{^f)-u{xt)) (^^^^^-tH^ff^^,.^^ ^^2) 

Using such a connection, it is immediate to obtain an expression of the conditional probability ([12]) in the form of a Feynman 
path integral 



P{xf,t\x,) = .-|(^(-/)-^(-.)) ^ r^'^=^^ 2,^ ,-/o<'^ ( &+v^fM )^ 

J x{ti)=Xi 



(13) 



where fA(^ is a normalization factor which comes from the Wiener measure and assures that / dx P{x^t\xi) = 1. The expression 
([13]) could have been obtained directly by computing the probability of path generated by iterating a discretized representation 
of the Langevin Eq. ^ — see e.g. the discussion in | 5, 15] — . The advantage of the derivation given here is that it does not 
involve the stochastic calculus. 

Eq. ( [T3] ) provides a microscopic representation of the conditional probabilities, formulated in terms of the Langevin trajecto- 
ries in configuration space which connect Xf and Xf.ln the next sections we shall discuss how such conditional probabilities can 
be effectively evaluated using the DRP formalism. 

The path integral expression of the conditional probability current ^ is 

J{xf,t\xi) = ^ {{v{x))[^ -D pV[/(x)) P{x,t\xi), (14) 
where {v{x))[. denotes the average velocity the system reaches the configuration Xi at time t, i.e. 

, f^^xx(t)e-^^ff^^^ 
Mx))l = ^-^ ^ , . (15) 
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III. THE DRP FORMALISM 



The DRP approach is based on the saddle-point approximation of the path integral ([13]). The expansion parameter controlling 
the accuracy of such an approximation is the thermal energy ksT , which enters in the definition of the diffusion constant D and 
of the effective potential y^.//(x). 

The saddle-points are the paths with the highest statistical weight Qy.^{—Seff[x]), i.e. those which minimize the effective 
action functional 

Seff[x]= f^dx(^^+Veff[x]Y (16) 

Note that all such so-called dominant reaction pathways satisfy the boundary conditions 

x(0) = Xi 

x{t) = Xf. (17) 
By imposing the extremum condition SSeff[x] = 0, we obtain the equation of motion for the dominant reaction pathways x{t)\ 

^^t) = VVeff[x{t)]- (18) 



In principle, a solution of the boundary- value problem ( 17)-( 18) may be obtained by minimizing numerically a discretized 



representation of the effective action functional Seff[x]. In practice, however, the presence of decoupling of time scales makes 
such a task very challenging. Indeed, computing a single dominant pathway would require to find a minimum of a function with 
a large number of degrees of freedom, dxNt, where Nt the number of time discretization steps. 

Fortunately, a major numerical simplification of this problem can be achieved by exploiting the fact that the equation of motion 
([18]) is simplectic, i.e. it conserves the "effective energy" 

Eeff=^iHt)-Veffm]- (19) 

Hence, rather than minimizing directly the effective action Seff[x], it is possible to obtain the dominant reaction pathways using 
the Hamilton- Jacobi (HJ) formulation of classical mechanics. In other words, the trajectories obeying the equation of motion 
( [T8] ) and subject to the boundary conditions ( [TT] ) are those which minimize the effective HJ functional 

ShMI)] = j'J dl^Eeff{t)^Veff[x{l)l (20) 

where dl = V dx^ is the measure of the distance covered by the system in configuration space, during the transition. The 
advantage of the HJ formulation is that it allows to remove the time as an independent variable. Instead, one introduces the 
curvilinear abscissa /. Since there is no gap in the length scales of molecular systems, the discretization of the HJ is expected to 
converge extremely much faster than the time discretization of the effective action Seff[x]. 

The effective energy E^/f is an external parameter which determines the time at which each configuration of a dominant 
reaction pathway is visited, according to the usual HJ relationship 



t{x)= r dl ' _ ^_^^ =. (21) 



1 

\/4D{Eefi+Veff[x{l)]y 



Typically, one is interested in studying transitions which terminate close the local minima of the potential energy U (x) . The 
residence time in such end-point configurations must be much longer than that in the configurations visited during the transition. 



From Eq. (21 ) it follows that these conditions are verified if 

Eeff^-Veff{Xo). (22) 

where Xo is a configuration in the vicinity of a local minimum of U {x) . 

In practice, computing the dominant reaction pathway connecting two given configurations Xf toxf amounts to minimizing a 
discretized version of the effective HJ functional: 



Ns-l 



(0] = L ^^[Eeff^yeff{x{n))] Mn^n^,, (23) 
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where Ns is the number of path discretization sHces. Once such a path has been determined, one can reconstruct the time at 



which each of the configurations is visited during the transition, using Eq. (21 ). In particular, the time interval between the ^-th 
and the {n + l)-th slice is 

A?„+i „ = A/„+i,„ 

^4D{Eeff + Veff[x{n)]y 

where A/^+i^^ = {x{n-\- 1) —x{n))'^. For a discussion on how to efficiently perform the relaxation of the HJ action, we refer 
the reader to liTOl[TTl, where the DRP method is used to investigate ab-initio chemical and conformational transitions in realistic 
molecular systems. 



Note that, in general, the solution of the boundary value problem ([T8])-([T7]) is not unique. Hence, one should in principle 
take into account for the entire set of dominant paths Xi{l) which obey the same boundary conditions ( 17 ). In practice, in many 
transitions of interest the relative statistical weight of secondary dominant paths with the same boundary conditions is much 
smaller and can be neglected. 

The dominant paths which solve the saddle-point equation ([18]) can be used to estimate the time evolution of an arbitrary 
configuration-dependent observable 0{x), during a transition from the reactant to the product. To this end, let hR{x) and hp{x) 
be the characteristic functions of the reactant and product states respectively — i.e. = 1 if x G R{P) and hi^(^p){x) = 

otherwise — and let po(-^) be the initial distribution of configurations in the reactant state. The average value of the observable 
0{x) at some intermediate time <x <t, evaluated over all possible reactive pathways which visit the product state at time t 
reads: 

fdxfhp{xf) !dxihR{xi) e-f(f/(-/)-^(-.)) po(^.) 2)xO[x(T)] e-V/W 

{0(t)) = 3 ', ,_ . (25) 

fdxfhp{xf) JdxihRixi) e-m^^/y^M) poixi) QZI^ ©xe-V/M 

To lowest-order in the DRP saddle-point approximation (i.e. up to corrections of order ksT) this average can be approximated 
with an average along the dominant paths only. In such an approximation, using the relationship Seff[x] = —Eefft-\-SHj[x], we 
have 

^ fdxf hpjxf) Jdxt hpjxt) ,-|(^(^/)-^(^.)) po(x,) I, 0[x,{x)] /^V^-^^^M ^^^^ 
Jdxf hp{xf) JdXi hR{xi) po(^.) e^'eff'-SHjM 

where the sum ^^^^ ^^^^ the dominant paths Xk{t), with boundary condition Xk{G) = Xf, Xk{t) = Xf. In principle, the 
effective energy parameters E^j^j^ have to be fixed in such a way that the total time of the transition is the same for all pathways. 
In practice, the effects arising from choosing the same value of Ee/f for all reaction pathways are usually found to be negligibly 
small. 

Once a dominant path x{x) has been determined, it is also possible to identify the configuration xjs along this path which 
belongs to the transition state (TS). This can be defined as the set of all configurations from which the system diffuses with 
probability 1/2 into the product, before visiting the reactant. Clearly, the reactive pathways cross the transition state. In par- 
ticular, the configurations of the dominant reaction paths which are representative of the TS can be found by requiring that the 
probability to diffuse back to the initial configuration Xi along the saddle-point path, equates that of evolving toward the final 
configuration. To the leading-order in the saddle-point approximation, this condition leads to the simple equation [8] : 



U{xf)-U{xi) 



[ dl ^J^{Eeff^Veff[xm - dl ^J^{Eeff^Veff[xm' (27) 



This equation can be easily solved for xjs^ once the dominant path x(x) has been calculated. 



IV. ACCOUNTING FOR SMALL FLUCTUATIONS AROUND THE DOMINANT PATHS 



Let us now go beyond the lowest-order saddle-point approximation, and compute the contribution of the stochastic fluctu- 
ations around the dominant paths, to order UbT . For sake of simplicity, in the following we shall consider the case in which 
the path integral associated to the conditional probability P{xf^t\xi) has a single saddle-point. The generalization to multiple 
dominant reaction pathways is straightforward: one simply needs to repeat such a calculation for each local maximum of the 
path probability. 

Let us consider a completely general transition pathway x(x) with boundary condition x{t) =Xf and x{0) = Xf and re- write it 
as the sum of the dominant trajectory x{z) and a fluctuation 3;(x) around it: 



x{x)=x{x)^y{x). 



(28) 
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Note that, by construction, the fluctuation path 3;(x) satisfies the boundary conditions y{0) =y{t) =0. 

We now provide an approximation to the path integral defining the conditional probability P{xf^t\xi) hy functionally expanding 
the action around x(x) : 

^ Seffix] + \jldx 1^ d^yiix) [x] y,(t), 



(29) 



where we have introduced the so-called fluctuation operator defined as 



?>Xi{x')?>Xk{x") 



1 (p- 

lb ^il^-^+^i^kVeffixit) 



5(x-x'). 



(30) 



Note that the fluctuation operator F[x] determines the amplitude of the stochastic fluctuations around the dominant path x{x). 
We also stress that the expansion of the effective action does not contain linear terms in the fluctuation field ^(t), because the 
dominant path x{x) around which we expand is a solution of the equations of motion ( 18 (, i.e. a stationary point of the effective 
action functional, Seff[x]. 

The functional integral over the fluctuation field y{x) can be performed formally. To this end, expand the fluctuation function 
y{x) in the basis of the (real) complete set of eigenfunctions of F[x]: 

y{'^)=Y, ' ^ I-^] = ^« ^" ■ ^^^^ 

n 

where the eigenfunctions x„ are chosen so as to satisfy the boundary conditions x„ (0) = x„ (f ) = and the normalization condition 



JO 



To second order in the fluctuations, the action reads 



^eff[- 



^eff[- 



dCn 



and the measure of the functional integral can be re- written as 
Hence, the conditional probability can be written as 

p{xf,t\xi) ^ r^^"' 

Jx{0)=Xi 
-\(U(xf)-U(xi)) 

= g\[ - ^ ^-V/W. 

detF[jc] 



(32) 



(33) 



(34) 



(35) 



(36) 



Eq. ([35| does not yet provide a practical tool to compute the conditional probability, since both the normalization factor 9^ 
and the determinant of the fluctuation operator detF[x] diverge in the continuum limit. However, there exist a standard trick to 



represent the ratio 



/detF 



in a form which is finite and calculable in the continuum limit. The idea is to introduce a so-called 



regulator, i.e. to multiply and divide by the conditional probability density of a fictitious system Freg.{xf-,t\xi). Such as system 
must be chosen in such a way that (i) the solutions of the Smoluchowski equation are known analytically and (ii) the DRP 
expression ( [35] ) to second-order in the fluctuations provides the complete exact result. 

For example, one can use as regulator the conditional probability associated to the diffusion in an external harmonic potential 
such as 



Uho{x) 



\ 



a {x-xo)^. 



(37) 
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The solution of the Smoluchowski equation for such a simple system are known analytically. In particular, in appendix [A| we 
show that for X/ = Xi = xq and choosing the parameter a in such a way that a^Dt ^ 1 one has: 



Preg.(xQ,t\xo) 



271 / 



d/i 



(38) 



It is immediate to verify that, for such a system, all the contributions to the expansion ( [29] ) beyond the second order in the 
fluctuation field j(x) vanish identically. Hence, the second-order DRP approximation yields in fact the correct exact result. 
Alternatively, one may use as regulator the conditional probability associated to the free Brownian motion: 



Preg.{xo,t\xo) 



') 



d/2 



4%Dt 



(39) 



Also for such a system, there is no contribution to the DRP expansion beyond the second order. 

In order to remove the divergences in Eq. ( [35] ), we multiply and divide by Preg. {^o^tl^o) and we replace the term in the 
denominator with its (exact) second-order DRP saddle-point representation: 

-l{U{xf)-U{xi)) 
PDRp{Xf,t\Xi) = . ^-'-//M 



detF[jc] 



Preg.Mxo) 



PregX^O^t\xo) 



Preg.{xQ,t\xQ) 



/detF[x] /T-^^ 



dQiF[x] 



det {Freg.[Xreg] 



exp 



hu{xf) - U{xi)) ^Eefft[x]-SHj[x] - \ Trlog (^-/.[x-,,,.] F[x]) 



(40) 



where E^Jr is the value of the effective energy parameter for which the total time in the regulator conditional probability is the 



^eff 

same as in the probability PDRp{xf^t\xi) we want to compute, i.e. 

r^/ di m 

Ix, 



t[x] = r 

J Xi 



^AD{Eeff + Veff[x\, 



dl 



^ '^^e/f i^reg.]) 



'-reg. [-^reg.i 



(41) 



eff 



Some comment on the expression (40) are in order. First of all, we observe that the factor det {Pr^l [xreg] F[x]) remains finite 
in the continuum limit. Then we note mat the o(kBT) correction in the DRP formalism is the analog of the Ginzburg correction 
of statistical field theory, and of the one-loop correction of quantum field theory. Finally, we emphasize again the fact that the 
regulator does not need to have any physical interpretation. It has been introduced as a mere mathematical trick to regulate the 
divergences appearing in Eq. ([35 ). 

From Eq. ( [4Q| ) it is straightforward to obtain the DRP expression for the probability current. 



J{xf,t\xi) = -D iy ^^'S/U{x))P{x,t\xi 
= D 



D 



(42) 



where uq {x) is the versor tangent to the dominant path at the configuration x. 

We note that in deriving Eq. ( [42| ), we have neglected the contribution coming form the gradient of the fluctuation determinant, 
since these term provides corrections which are of higher order in ^^T. The o{kBT) DRP expression for the probability current 
will be used in section [Vll| to compute the reaction rates. 

The DRP formula for the time evolution of the average of the observable 0{x) including o{kBT) corrections reads: 



fdxf hp{xf) Jdxi hR{xi) e-m^f)-"(^^)) po(x,) I, 0[x,(x)] /^//"^"^[^"d /^et (F,,g.[xo] F-i[x,]) 

{0{i)) ^ 1 7 . ^ . (43) 

!dxf hp{xf) fdxi hR{xi) po(^.) /.//-^"^[-dJdet [xo] 
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FIG. 2: The total time for a transition from Xi = —I to Xf 
discretization steps, in the original theory and in the EST. 



1 in the quartic potential pSl), calculated using Eq. (I24l with different number of 



As in Eq. (26 ), the sum Y,k ^^^^ over all the dominant paths Xk{t). 



In practice, the determinant det {Freg.[xreg.] F ^[x]) has to be evaluated numerically from a discretized representation of the 
fluctuation operator associated to the dominant path F[x]. To obtain such a representation, one needs to express the time deriva- 
tive using discretized time intervals. It is most convenient to use the intervals Ati^i^i evaluated from the dominant trajectory, 
according to Eq. ([24]). The fluctuation operator reads 



Fi^yJ - zllR 



At 



m+l,m 



At, 



m+l,m 



At, 



m,m— 1 



At 



m,m— 1 



(44) 



where the indexes m - 
system. 



^Ns run over the path frames, while the indexes ij= 1, . . . , label the degrees of freedom of the 



V. IMPROVING THE CONVERGENCE ON THE DRP CALCULATION 



The main numerical advantage of the DRP formalism arises from the possibility of replacing the time as a dynamical variable, 
and replace it with the curvilinear abscissa /. Since there is no gap in the length scales of molecular systems, the convergence 
of the discretization of / is usually very fast. As a result, using such a formulation, it is possible to gain information about the 
reaction mechanism at a very low computational cost. 

On the other hand, in order to obtain information about the dynamics, one needs to compute the times at which each configu- 
ration is visited along the dominant path, using Eq.s ( [2T] ) and (24 ). The calculation of the time intervals A^^ ^+1 is also needed in 
the discretized representation of the fluctuation operator ( [44] ), which enters in the calculation of the order ksT corrections arising 
from non-equilibrium stochastic fluctuations around the dominant path. 

A potential limitation of the DRP approach resides in the fact that the DRP calculation of the time intervals converges very 
slowly with the number of equal displacement discretization steps. As a result, in order to achieve an accurate description of the 
dynamics, or in order properly take into account of the effects of fluctuations, one needs to use a large number of path frames, 
with a consequent significant increase of the computational cost of a DRP simulation. 

To illustrate this problem in a simple example, let us consider the diffusion of a particle in a one-dimensional quartic external 
potential 



2\2 



U{x) a{x{)-x) 



(45) 



We shall adopts units in which = \, ^ = \ / ksT = 5 and Z) = 1. 

The dashed line in the left panel of Fig [2] shows the result of the calculation of the total time interval for a dominant transition 
from X; = —1 to X/ = 1, as a function of the inverse of the number of equally-displaced path discretization steps Ns. Such a 
time interval was computed by adding up all the elementary time intervals A^^ evaluated according to Eq. (24), by choosing 
Eeff = —Veff{xi) +0.03. This figure shows that, in order to reduce the discretization errors below 1%, one needs to use ~ 10^ 
path discretization steps. 
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Such a slow convergence is of course a consequence of the decoupling of the time scales characterizing the dynamics of this 
system. Indeed, the quasi-free diffusion in the bottom of the wells is much slower than the crossing of the transition regions, 
where the force is large. An analog convergence problem is encountered also in molecular systems, since the characteristic 
internal time scales are decoupled and range from fractions of ps to ns. 

The purpose of this section is to show that the convergence of the calculation of the time intervals in the DRP approach can 
be greatly improved by adopting the effective stochastic theory (EST) developed in llT4l and briefly reviewed in the appendix |B] 

The main idea of the EST is to exploit the gap in the internal time scales in order to analytically perform the integral over the 



fast Fourier components of the paths x(x) which contribute to the path integral ( 13 ). Through such a procedure, the effects of the 



fast dynamics is rigorously and systematically averaged out by "renormalizing" the effective potential: 

Veff{x) ^ V^ff{x)=Veff{x)+V^ff{x). (46) 

^///W = ^^|ilr^vV.,,(x) + ... (47) 

In such an Eq., Ar^ is a cut-off time scale which must be chosen much smaller than the fastest internal dynamical time scale and 
is a parameter which defines the interval of Fourier modes which are being analytically integrated out — see the discussion in 
the appendix |b] — . Typically, for molecular systems A^^ ^ 10~^ps and b ^ 10~^ ifTSTl . The dots in Eq. (47 ) denote higher order 
correction in an expansion in the ratio of slow and fast time scales (slow-mode perturbation theory). 

The EST generates by construction the same long-time dynamics of the original — or so-called "bare" — theory, but has 
a lower time resolution. In the context of MD simulations, this implies that the EST can be integrated using much larger 
discretization time steps L15 ] . In the context of the DRP simulations, the utility of the EST resides in the fact that fewer path 
discretization time steps are required in order to achieve a convergent calculation of the time interval from the dominant path, 
through Eq. (21 ). Such a gain is clearly visible in Fig. [2j where we compare the total time interval obtained in the bare theory 



and in the EST, for different numbers of path discretization steps Ns. From the mathematical point of view, the computational 

^EST 
^eff 



gain of adopting the EST can be seen as a consequence of the fact that the renormalized effective potential V^fl {x) is in general 



a smoother function than the bare effective potential Veff{x). 



VI. TESTING THE DRP CALCULATION WITH o{kBT) CORRECTIONS ON AN ANALYTICALLY SOLVABLE MODEL 

In this section, we assess the accuracy of the DRP o{kBT) calculation developed in the previous section by computing the 
conditional probability and the probability current in an exactly solvable model. In particular, let us consider a one-dimensional 
point particle diffusing in a harmonic oscillator of potential 

UHO = loix\ (48) 

The conditional probability for the point particle to be at the origin xi = at the initial time and to reach the point Xf after a 
time interval t is known analytically and reads: 



M oB f OL^xl cosh(a ^Dt) I ] 

PHo{xf,t\xi = 0) = J-- . , , V> ^ , exp<^ , { , , ' ^ , ^ -a^Dt} (49) 

^ ^' ' ^ ]l 4% smh{a ^ D t) ^1 4smh{a^Dt) 2 ^ J 

Let us now discuss the DRP calculation of the same conditional probability. We choose to use as regulator the conditional 
probability of the free Browian diffusion for Xf = Xf, i.e. 



^o(xnfk) = ;/^ (50) 



The NLO DRP result is therefore 



where Fho [x] is the fluctuation operator of the harmonic oscillator evaluated along the dominant path x{x) evaluated numerically, 
Fo[xi] is the fluctuation operator of the free Brownian theory, evaluated on the static path x{t) = Xf — with Seff[xi] = — , while 
t[x] and Shj[x] are computed from the dominant path x{x) using Eq.s (20) and (21 ), respectively. 

We recall that, in the specific case of the diffusion in an harmonic oscillator, all the contributions to the DRP saddle-point 
expansion beyond the second order vanish identically. Hence, the NLO prediction ([5T]) must agree to numerical accuracy with 
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FIG. 3: Comparison between the DRP expression for the conditional probabiUty for the diffusion in a harmonic potential, in a system of units in 
which a = 4 P = 10 and = 1 . In the left panel we show the DRP and the corresponding exact conditional probability P{xf,t\xi) as function 
of the final position Xf, for Xi = and for two times t = 0.1 and t = I. We also show the corresponding Boltzmann equilibrium distribution. In 
the right panel, we compare the DRP prediction for the probability density and for the probability current with the corresponding exact results, 
as a function of the time interval t, at a fixed final position Xf = 0.2. 



the analytic result ( [49| ), for all choices of t and Xf. In the left panel of Fig. plwe compare the DRP and the exact analytic 
conditional probability for Xi = as function of the final position Xf, at the fixed times ^ = 0.1 and ^ = 1 for P = 10. We note 
that at long times — i.e. t = I — the system has attained thermal equilibrium (i.e. the Boltzmann distribution). In addition, 
the DRP calculation yields the exact result at any intermediate time, as expected. In the right panel of Fig. [3] we compare the 
DRP prediction for the probability density and for the probability current with the corresponding exact results, as a function 
of the time interval t, at a fixed final position Xf = 0.2 for P = 10. We see that, in the long time limit, the probability current 
progressively dies out as the corresponding probability distribution becomes stationary. As in the previous case, the result of the 
numerical DRP calculation completely agrees with the expected analytic prediction. 



VII. REACTION KINETICS 



In the previous sections we have discussed and tested the numerical evaluation of the conditional probability P{xf^t\xi) 
and of the probability current J{xf,t\xi), in the DRP approach. In this section, we show how such quantities can be used to 
microscopically compute the reaction rate. 

Let us consider the case of a conformational reaction involving two thermodynamically (meta)-stable states, which we shall 
refer to as to the reactant R and product state P, respectively. The thermodynamical states are usually defined in terms of an 
order parameter O, i.e. a quantity which is distributed around different values in the reactant state R and in the product state P, 
i.e. O Or in R and O Op in P. Hence, along the reaction pathways, O varies very rapidly from Or to Op. 

We shall further assume that the time scales in which the thermalization is achieved in each of the two-states is much smaller 
than the inverse rate of transitions between them. Under such conditions, the system is said to obey two-state kinetics. In this 
case, the fraction of population in the reactant and product obey the well-known kinetic Eq.s: 

nR{t) = -kp^p nR{t) ^ kp^R np{t) (52) 
rip{t) = kp^p nR{t)- kp^Rnp{t). (53) 

with the condition np{t) -\-np{t) = 1. 

At equilibrium, the population fraction stops depending on time, (t) = n^^p^ , and Eq.s ( 52 )-( 53 ) yield the detailed balance 
condition, 

^ = (54) 

rip kp^p 



The solution of the kinetic Eq.s (52)-(53) with boundary conditions np{0) = 0^np{0) = 1 — i.e. assuming that the system is 
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initially prepared in the reactant state — are 

np{t) = n^{\-e-^'^, (55) 
nR(t) = l-np{t), (56) 

where k = kp^p + kp^p is the thermal relaxation rate. 

In the opposite short-time regime, i.e. for ^ <C p np{t) grows linearly with the reaction rate kp^p: 

np{t)c^kR^pt {kR^pt<^l). (57) 

Hence, the reaction rate kp^p can be written as 

kp^P = lim^np{t). (58) 
t^o at 

The limit in such an equation means that the time t is chosen much smaller than the inverse reaction rate — i.e. kp^p ^ <C 1 — , 
yet much larger than all the thermalization time scales in the reactant state. The possibility of making such a choice is guaranteed 
by the assumption of two-state kinetics. 

In order to microscopically compute the rate in the DRP approach, we need the path integral expression of the population 
fraction np{t), solution of the Eq.s ( 52 )-( 53 ) with boundary condition np{0) = 0. In the appendix[c]we show that such a solution 
is given by 



np{t) = J dxf hp{xf) J dxi hp{xi) po(x,-) P{xf,t\xi) 



91 fdxfhp{xf) /jx,-/zHxOPo(^0^"^^''^"^^"''^"'-^^ r^'^ "^^^,-/o^<l^+%,W 

J J Jx{ti)=Xi 



(59) 



where po(-^) is the initial normalized distribution of the configurations in the reactant state and hp{x) and hp{x) are the charac- 
teristic functions of the reactant and product state, respectively. 
Eq. ( [58] ) becomes 

kp^p ^^^^=^^^ j dxfhp{xf) J dxi hp{xi)po{xi) ^^P{xf,t\xi). (60) 
Using the Smoluchowski Eq. we can rewrite this as 

kp^p ^^^=^^ Yt^^^^^^ I ^^f^P^^f^ j dxihp{xi)po{xi)^^P{xf,t\xi) (61) 

= J dxfhp{xf) j dxihp{xi)po{xi)V -J{xf,t\xi) (62) 

We now introduce a closed surface dW which surrounds the product state. Using Gauss's divergence theorem, we find: 

kp^p ^^^=^^ 1 dxi hp{xi) po{xi) 1 da-J{xf,t\xi) (63) 
J Jdw 



■I 



a| hx-J{xf,t\xi) (64) 

dw 



where J{xf^t\xi) denotes the average of the current with respect to the initial configuration in the reactant, i.e. 

J{xf,t\xi) = J dxi hp{xi) poixi) J{xf,t\xi), (65) 

and hx is the unitary vector normal to the dividing surface dW at the point x, oriented out- ward. Note that if the surface dW is 
chosen in such a way that it intersects the transition region, then Eq. ([64]) yields in fact a multi-dimensional generalization of 
Kramers' flux-over-population expression for the reaction rate| 16 1. 

The Eq.([63]) contains the integrals over two large dimensional spaces, which may be difficult to perform in practical applica- 
tions. Hence, it is useful to introduce further approximations. First of all, we observe that under the assumption of two-state 
kinetics, the current J{xf^t\xi) becomes quasi-instantaneously independent on the initial configuration Xi. This fact allows to 
remove the average over the initial configurations in the reactant, since for any choice of Xi G R one has 



J(xf,t\xi) J{xf,t\xi) yxi e R. (66) 



12 




FIG. 4: The definition of dW as the hyper-surface orthogonal to the dominant path the the point xts solution of the transition state Eq. ^2j\ . 



Let us consider first the case in which the reaction occurs through a single channel, i.e. that all the stochastic paths starting form 
Xi and ending up in the product after a short time t are confined in a small bundle around a single dominant reaction pathway. In 
order to estimate the flux of the current through the dividing surface, we insert the DRP expression for the current, i.e. 

W ^ - / d\a\,h,-JDRp{x,t\Xi)^PregXxo,t\xo)e-^e^^^^^^^ 

Jdw 



In the low-temperature limit, the leading dependence on x in the integrand of Eq. ( 67 ) comes from the exponential factors and 
we can neglect the dependence on x of all the non-exponential terms. In addition, in the same limit, we can use the approximation 



Shj[^ 



l^£j/|Vt/[x-(/)]| = ^(t/(x)-[/(xO) 



2kBT 



Hence, up to corrections of higher order in ^^T, the flux ( [67] ) takes the form 

kR^P = - / d\a\j h^'J{x,t\xi) -Z^w e^^^'^Ts) ^fi • Jj)j^p{x^^,t\xi)), 
Jxedw 

where x^^r is a configuration on the dividing surface dW and Z^^r is the partition function 

Jxedw 



(68) 



(69) 



(70) 



Eq. ( |69| ) is independent on the specific choice of the dividing surface. We now specialize to the case in which dW is the 
hyper-surface orthogonal to the dominant reaction pathway, at the point xjs solution of the transition state Eq. ([27]) — see 
Fig. [4] — . With such a choice. 



mTs)[ 



(71) 



Note that uq{xts) is the unit vector tangent to the dominant path at the configuration xts (hs is the value of the curvilinear 
abscissa a for which Eq. ([27]) is satisfied). Correspondingly, the partition function reads 



Z^w=Zts = J dx5[{x-XTs)-UQ{^Ts)] e 



Pt/(x) 



(72) 



The partition function (72) can be estimated in local harmonic approximation, by running short MD simulations starting from 
xts^ subject to the constraint to lie on the surface orthogonal to the tangent to the reaction pathway at xjs^ i-e. 



Zjs - 



(73) 



13 



where x' is the /— th coordinate of the configuration x and 



{[X -Xjs) U = or/ X ^ 1 _Rr/rr^ • ^/^^^ 



Hence, using the DRP expression for the reaction current we arrive to our final result: 



^ weW - ^^f-fi^Ts) 



V^det(F..,.[x..,]f-l[x]) 



(75) 

In such an Eq., the effective energy parameter Ee/f must be chosen in such a way that the total time t[x] evaluated according 



to Eq. (21 ) is much larger than the relaxation time in the reactant and yet much smaller than the total relaxation rate. In such a 



time regime, the flux of probability current across the transition state is stationary and the value of the expression ( 75 ) must be 
independent on the specific value of Ee/f chosen. 

Finally, if the reaction can occur through more than one reaction pathway, the total rate is obtained simply by adding up all 
such contributions: 

W ^ XZ^,.P^(45) \J{x'rs.t\xi)\. (76) 

k 

For very simple systems, computing the rate by means of Eq.([75]) of Eq. ( [76| ) is expected to be more computationally expensive 
than in standard Kramers theory 1 16 |. On the other hand, the advantage of the DRP method developed here is that it does not 
require to know a priori the location of the transition state. Hence, we expect that our method may be used to investigate the 
kinetics of two-state reactions in large configuration spaces, which are generally characterized by a complicated energy surface. 
In particular, for many complex molecular systems, the transition state cannot by guessed from the structure of the interaction, 
and the multi-dimensional Kramers theory is therefore useless. 

We also observe that the DRP method presented in this section bares some similarity with other existing techniques for rate 
calculation. In particular, the identification of the reaction coordinate / from a statistical important reaction pathway is also 
used in the so-called milestoning method 1 17 |. An important advantage of such an approach with respect to the DRP method 
developed here is that it does not require to assume two-state kinetics. On the other hand, the present method is much less 
computationally expensive, as it does not require to evaluate the first-passage-time distributions from the different milestones, 
from MD simulations. 



The DRP Eq. ( [75) bares also some similarity also with Chandler's theory 1 18] and with transition state theory 1 19 |. Indeed, 
in both such approaches, the rate is related to the flux of reactive trajectories through the transition state. On the other hand, we 
stress the fact that in the present DRP approach, such a flux is evaluated over non-equilibrium trajectories. 

Finally we note that, in the transition path sampling algorithm |2|, the rate is usually calculated starting from the path integral 



expression of the product population fraction, i.e. Eq. ( [59] ). However, in such an approach, the normalization of the path integral 
is obtained by evaluating the free energy of the transition path ensemble, i.e. the reversible work which is required to constraint 
the final configurations of the paths into the product state. On the other hand, in the DRP approach, such a normalization is 
guaranteed by construction a priori, by the regularization procedure. 

VIII. TESTING THE DRP CALCULATION OF THE REACTION RATES 

In order to test the scheme developed in the previous section for rate calculations, we study the reaction kinetics of some 
simple two-dimensional toy systems, for which very accurate results can be obtained also using other methods. In particular, in 
wide range of temperatures, the rate for such systems can be accurately evaluated using Kramers theory or computed directly by 
running long MD simulations. 

Let us begin by considering the rate of escape from a standard two-dimensional bi- stable potential, 

U{x,y) = a{x^ - ay^)^^by^ (11) 

with a = 1/8, b = 2, and CO = 2. In this simple case, the dominant path connecting = — CO to x/ = CO is known analytically 
(straight horizontal line). We have represented such a path using Ns = 150 equally spaced steps. In addition, performing an 
accurate calculation of the fluctuation determinants for such a simple two-dimensional system is straightforward, even without 
residing on the low time resolution effective description. 

We recall that in the DRP approach, the total time of the transition is determined by the value of the effective energy parameter 



Eeff, according to Eq. (21 ). The same parameter is used in the calculation of all the discretized time intervals which enter in the 
fluctuation determinant. On the other hand, our prediction for the rate must obviously not depend on the choice of Ee/f. 
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FIG. 5: Left panel: The absolute value of the normalized probability current at the sadde-point Xf = 0, evaluated with the DRP method at 
different times for = — CO and P = 5. The different total times are determined by different choices of the effective energy parameter. Righ 
panel: The ratio between the reaction rates evaluated with the DRP method and with Kramers theory. 



In order to see how this condition can be satisfied, we recall that our rate expression (58) requires the time interval t to be 
much smaller than the inverse rate, kR^pt <C 1 . At the same time, t must be chosen much larger than the thermalization time 
in the reactant state, to assure single exponential relaxation. If both conditions are simultaneously satisfied, then one should 



observe a stationary flux through the transition state. Consequently, our rate expression (76) should become independent on the 
precise choice of ^ — hence of Ee/f — . 

In order to test if such a stationary current is realized in our simulations, in the left panel of Fig. [5] we plot the DRP current at 
the saddle-point as a function of the total time (hence for different values of the effective energy). We can clearly see that for the 
longest times the current becomes stationary, hence the rate stops depending on the specific choice of Ee/f. The right panel of 
Fig. [5] shows that, in this system, the calculation of the rate using Kramers' theory using the DRP approach agree within 1% 
accuracy. 

Let us now consider a toy model which allows to specifically assess the role played by the stochastic fluctuations around the 
dominant path, in the kinetics of the reaction. To this end, we study the two-dimensional three-state system, consisting of a 
reactant R, and of two products Pi and P2, shown in the left panel of Fig. [6]The functional form of the potential energy is 



U[x,y] :-- 



-bx" 



- k\ exp 



2a 



+ (/-2x^ + ia: + ^l)/ 



(78) 



with a = 5/64, b = —10/16, c = 1, ^1 = —0.6, xi = 0, a = 0.3, k = 0.5 and Q. = 2. Also in this case, the dominant paths 
are known analytically. The HJ action and the fluctuation determinants appearing in Eq. ( [75] ) are evaluated using A(y = 100 
discretization steps. 

This model is built is such a way that the energy barrier separating the reactant to the two products is the same. Yet, the steeper 
structure of the potential energy on the right of the reactant tends to disfavor large stochastic fluctuations, along the dominant 
path of the reaction R ^ P2. We therefore expect that the rate of escape from R to Pi should be larger than that from Rio P2. 
This fact is clearly seen in the right panel of Fig. [6j were we compare the ratio kp^p^kR^p^ evaluated using Kramers theory 
and using the DRP method. We see that the two approaches give results which agree within about 2% accuracy. In both cases, 
the reaction P ^ is about 50% slower than the reaction P ^ P\. 



IX. CONCLUSIONS 



In this paper, we have analyzed the role of the stochastic fluctuations around the most probable reaction pathways, in systems 
obeying the over-damped Langevin dynamics. Such fluctuations affect the probability for a transition to take place through a 
specific reaction channel in a given time, hence the kinetic of the reaction. For sufficiently small temperatures, the fluctuations 
in the ensemble of reaction pathways are confined within small bundles around the locally most probable reaction pathways. In 
such a regime, we have developed a technique to efficiently compute their contribution, to order ksT accuracy. Clearly, in the 
opposite high-temperature regime, the stochastic fluctuations become very large and overshadow the information encoded in the 
dominant paths. In this case, the accuracy of the DRP approach breaks down. 
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We have shown that the calculation of the time intervals and of the fluctuation determinant can be made much more efficient 
by adopting a low-time resolution effective description, in which the fast dynamics is analytically pre-averaged out. Indeed, in 
the effective stochastic theory, the convergence of DRP calculations is achieved using much fewer discretization paths. 

The o{kBT) expression of the conditional probability was used to derive a formula for the reaction rate, within a flux-over- 
population approach. We have illustrated and tested our results on simple toy models, where accurate or even exact results could 
be obtained with alternative methods. In the future, we plan to apply these techniques to investigate the kinetics of much more 
complicated molecular reactions. 
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Appendix A: The Conditional Probability for the Diffusion in a J— Dimensional Harmonic Potential 



Here we review the exact calculation of the conditional probability for a point-particle diffusing in the J— dimensional har- 
monic external potential 



(Al) 



For sake of simplicity, we discuss such an evaluation in the case Xf =Xi = and we choose the origin of the configuration space 
in such a way Xq = 0. We also choose units in which the viscosity y = 1/pZ) is set to 1. 
In this case the effective action of the static solution x = is Seff[x] = — | CO?, hence: 

PHo{0,t\0) = CSCei ^> 



det 



P 32 , |3„4>, 



8o( 



(A2) 



det 



5^7(-£ + ®f) 



where in the last step we have absorbed the factor P/2 which appears in the fluctuation operators into the normalization constant 
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The normalized eigen-functions and eigen-values of the fluctuation operator with the boundary conditions y{t) = y{0) = are: 



sin 



n n 



(7ln)2 



-oof (n = 1 , 2, . . . , i = I, . . . ,d) 



Notice that the same eigenfunctions are also eigenstates of the fluctuation operator for the free diffusion 

32 



with eigenvalues 



Fo = -dijd{x'-x) 



(7) 1>2,..., /= l,...,d) 



8x2' 



(A3) 
(A4) 

(A5) 
(A6) 



As usual, in order to get rid of the unknown normalization 9{J we multiply and divide by a regulator, in this case the free 
propagator Po(0,f|0): 



PHo{0,t\0) = Po(0,r|0)e2 ^- 



nil n„>i ) 



471? 



0^ 



(A7) 



(A8) 



Using the result 



we find 



T-rf. I (i^ft \ ^ Sinh(o),0 



cor 



sinh(0L)?r) 



(A9) 



(AlO) 



Notice that, in the long time hmit, sinh(co/^) ^ so P//o(0,^|0) converges to the inverse of the partition function, as it 

should: 



;,„„(o,,io)^(|)'"'n<«,=^ 

Notice also that, for this system, the thermalization time does not depend on the temperature. 



(All) 



Appendix B: Effective Stochastic Theory 

In this section, we sketch the derivation of the EST. For all further details we refer the reader to the original paper lfT4ll . For 
simplicity and without loss of generality, it is convenient to consider the path integral with periodic boundary conditions 



Z{t) = J dxP{x\x\t) = 

The starting point to develop the EST consists in introducing the Fourier components of the paths, 

x((On) = - f dxx(x)e-''^-' 
t Jo 

x{x) = x{x^t)=Y,4(^n)e''''^'. 



(Bl) 

(B2) 
(B3) 
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where co„ ^ n, are the Fourier frequencies (fr = 0, ±1 , ±2, . . .) . 

The path integral ( |B1| ) is defined in the continuum Hmit. Numerical simulations are always performed using a finite discretiza- 
tion time step At. Clearly, the shortest time intervals which can be explored in a numerical simulation is of the order of few At. 
Equivalently, the largest frequencies of the Fourier transform of the stochastic paths x(co) are of the order few fractions of an 
ultra-violet (UV) cut-off Q = 2k /At. 

Let us now split the Fourier modes of the paths contributing to ([13]) in high-frequency — or "fast" — modes and low-frequency 
— or "slow" — modes. To this end, we introduce a real number < Z? < 1 such that the frequency range (0,^1) is split in two 
intervals (0, bQ,)[J {b ^1, Q,). Correspondingly, one can define the "fast" component of the path x> (x) and the "slow" component 
of the path x< (x), by summing over the Fourier modes in the (0, b Q.) and {b ^1, Q.) range, respectively: 

x<(0 = ^ x((0„)e'""' (B4) 

ba<\((in\<^ 



The complete path integral (Bl ) can therefore be exactly re- written in the following way: 
where 



is called the renormalized part of the effective action. 

The EST is constructed by explicitly evaluating Sy [x<], i.e. by performing the path integral over fast modes x> (x). In the limit 
in which the fast and slow modes are separated by a large gap in the spectrum of Fourier modes — i.e. if the system displays a 
decoupling of time scales — such an integral can be carried out analytically in a perturbative approach based on Feynman diagram 
techniques ifni . The expansion parameter such a perturbation theory is the ratio between the typical frequency CO of the slow 
modes and the UV cut-off bQ.. Clearly, if hard and slow modes are decoupled, the ratio co/ (bO.) is a small number, hence the 
terms proportional to higher and higher powers L of such a ratio provide smaller and smaller corrections. 

If one accounts only for the leading corrections in the l/bQ, expansion, the renormalized part of the action takes the form of 
an effective interaction term [14 1, i.e. 

where 



We emphasize that the result of the EST construction is a new expression for the same path integral (Bl ), in which the UV 
cutoff been lowered from to bQ.. Equivalently, the path integral is discretized according to a larger elementary time step. 
At At/b: 

Z^{t) ^ /(Dx^-V/Woc/ ^xe-'^ff^'^-^^'''''ffl<^^^ (BIO) 

J At J At/b 

In these expressions, the symbol denotes the fact that the path integral is discretized according to an elementary time step 
and we have suppressed the subscript "<", in the paths. It can be shown that the proportionality factor between Z^{t) and 
^Esri^) d^P^nds only on t and does not contribute to the statistical averages. 



Appendix C: Path Integral Expression for the Population Fractions np{t) and nR{t) 



In this section we provide a microscopic representation of the solutions np{t) and nR{t) of the kinetic Eq.s (52)-(53 ). Let us 
consider in particular the case in which the system is initially prepared in the reactant state, i.e. ^/?(0) = 1 and np(0)= 0. Let 
po{x) be the initial distribution of the configurations in the reactant state. Clearly, such a choice of initial conditions implies the 
normalization condition 



/ 



dx hR{x) po{x) = 1, 



(CD 
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where hR{x) is the characteristic function of the reactant, i.e. hR{x) = 1 if x e R and otherwise. Under the assumption of 
two state kinetics, the thermaUzation in the reactant occurs over a very short time scale, hence the specific choice of the initial 
distribution in the reactant is in fact irrelevant. 

We now show that a microscopic representation for the product population fraction np{t) is obtained by averaging the condi- 
tional probability P{xf,t\xi) over the initial configurations Xi in the reactant and summing over all possible final configurations 
Xf in the product, i.e. 

np{t) = J dxfhp{xf) J dxi hp{xi) po{xi) P{xf,t\xi) 

= fdxf hp{xf) f dxi hp{xi)M^i) r^'^^'' q)x ,-^o^<l^+%,w) ^^^^ 

J J J x{ti)=Xi 

We observe that Eq. (C2) satisfies the correct initial condition, ^p(O) = 0. 

Let us now introduce the complete set of eigenstates ^^(x) of the "quantum" Hamiltonian Heff, 

Heff^n{x)=kn'¥,{x). (C3) 

In particular, it is inmiediate to verify that the ground state of He// has a vanishing eigenvalue and reads 

-|t/(x) 

^oW = ^^, (C4) 

where Z is the partition function of the system, 

Z = J Jx^-P^W. (C5) 



By inserting the resolution of the identity, 1 =Y.n\^) (^1' i^t^ the "quantum" propagator ( 12) we obtain the so-called spectral 
representation of the conditional probability: 

P{x,t\xi) = ^-P/2(t/W-t/(x,)) £ ^^^{xi)e-''^'. (C6) 

Hence, the conditional probability P{xf^t\xi) converges to the Boltzmann distribution, in the long time limit: 

l.-P^W, (C7) 

regardless of the initial condition. 

In particular, the systems obeying two-state kinetics are those in which the spectrum displays a gap between the first and 
second eigenstates of the effective "quantum" Hamiltonian He/f'. 

h<^k2. (C8) 

Indeed, in this case the time scale = ^ decouples from all the other relaxation time scales in the system, and the approach to 
thermal equilibrium occurs through a single-exponential relaxation: 

P{x,t\Xi) +^-P/2(t/W-t/fe)) XI/t(^)xj;^(^.) ^-ht^ ^C9) 

If the reaction is two-state and if Xi and Xf are not in the same state (for example Xi G R and Xf e P) then the probability of 
performing a transition from Xi to Xf vanishes in the short- time limit, i.e. 

lim P{xf,t\xi)=0. (CIO) 

This fact implies that 

^5(x/)^o(^0 = -^I(x/)^i(xO, (Cll) 

thus Eq. ( |C9| ) gives 

P{xf,t\xi) = ^^ [l-e-'^^y (C12) 
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We emphasize that Eq. (CIO) — and therefore Eq. (Cll ) — are fT(9^ justified if Xf and Xf are in the same state. Indeed, in the 
approximation of two-state kinetics, the local thermalization in the R and P states is assumed to occur instantaneously, as the 
only finite time scale is the mean-first-passage time across the barrier. 
Using the spectral decomposition ( |C9| ) we find 

np{t) J dxfhp{xf) J dxihRixi) po(x,) 
The product population fraction at time t then reads 

np{t) = Hp ^1 



-k,t 



(C13) 



(C14) 



where Hence, we have recovered Eq. ( 55 ) and we have shown that first excited state of the quantum effective 

Hamiltonian is the equilibrium relaxation rate of the system, k\ = k. 
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